##########################################################
# John Henderson and Alex Theodoridis
# Replication Data for: "Seeing Spots", 
#  Forthcoming in Political Behavior, August 20, 2017
# 
##########################################################
#
#  figure_xii.R
#  -- file produces the mean bar plots in Figure xii
#
##########################################################

rm(list=ls())
source('~/Dropbox/Seeing_Spots/replication/pre_data_vagov.R')
library(gplots)
library(plotrix)


video_skipped=va_data$video_skipped_60             

ymat=cbind(video_skipped_60, replay, share, getlink, all_y_60,time_watched, total_time)   
   
is.even <- function(x) x %% 2 == 0         
# (a) pos v. neg
# outcome matrices for pos v. neg by pid   

outcomes=list()        
for(k in 1:7){
	outcomes[[k]]=matrix(NA,3,4)  
	colnames(outcomes[[k]])=c('Mean Positive','Mean Negative','Difference','pvalue')
	rownames(outcomes[[k]])=c('Democrats','Independents','Republicans')	
}
        
dpos=which(pid_lean==-1 & pos_ad==T); dneg=which(pid_lean==-1 & pos_ad==F)
ipos=which(pid_lean== 0 & pos_ad==T); ineg=which(pid_lean== 0 & pos_ad==F)
rpos=which(pid_lean== 1 & pos_ad==T); rneg=which(pid_lean== 1 & pos_ad==F)

for(k in 1:7){
	outcomes[[k]][,1] = c(mean(ymat[dpos,k],na.rm=T),mean(ymat[ipos,k],na.rm=T),mean(ymat[rpos,k],na.rm=T))
	outcomes[[k]][,2] = c(mean(ymat[dneg,k],na.rm=T),mean(ymat[ineg,k],na.rm=T),mean(ymat[rneg,k],na.rm=T))	
	outcomes[[k]][,3] = outcomes[[k]][,1] - outcomes[[k]][,2]                                              
	outcomes[[k]][,4] = c(t.test(ymat[dpos,k],ymat[dneg,k])$p.value,t.test(ymat[ipos,k],ymat[ineg,k])$p.value,t.test(ymat[rpos,k],ymat[rneg,k])$p.value)	
}

tone_by_pid=list('skipped'=outcomes[[1]],'replay'=outcomes[[2]],'share'=outcomes[[3]],'getlink'=outcomes[[4]],'all'=outcomes[[5]],'time watched'=outcomes[[6]],'total watched'=outcomes[[7]])
    
# (b) dem v. rep 
# outcome matrices for pos v. neg by pid   

outcomes=list()        
for(k in 1:7){
	outcomes[[k]]=matrix(NA,3,4)  
	colnames(outcomes[[k]])=c('Mean McAuliffe','Mean Cuccinelli','Difference','pvalue')
	rownames(outcomes[[k]])=c('Democrats','Independents','Republicans')	
}
        
dpos=which(pid_lean==-1 & dem_ad==T); dneg=which(pid_lean==-1 & dem_ad==F)
ipos=which(pid_lean== 0 & dem_ad==T); ineg=which(pid_lean== 0 & dem_ad==F)
rpos=which(pid_lean== 1 & dem_ad==T); rneg=which(pid_lean== 1 & dem_ad==F)

for(k in 1:7){
	outcomes[[k]][,1] = c(mean(ymat[dpos,k],na.rm=T),mean(ymat[ipos,k],na.rm=T),mean(ymat[rpos,k],na.rm=T))
	outcomes[[k]][,2] = c(mean(ymat[dneg,k],na.rm=T),mean(ymat[ineg,k],na.rm=T),mean(ymat[rneg,k],na.rm=T))	
	outcomes[[k]][,3] = outcomes[[k]][,1] - outcomes[[k]][,2]                                              
	outcomes[[k]][,4] = c(t.test(ymat[dpos,k],ymat[dneg,k])$p.value,t.test(ymat[ipos,k],ymat[ineg,k])$p.value,t.test(ymat[rpos,k],ymat[rneg,k])$p.value)	
}

cand_by_pid=list('skipped'=outcomes[[1]],'replay'=outcomes[[2]],'share'=outcomes[[3]],'getlink'=outcomes[[4]],'all'=outcomes[[5]],'time watched'=outcomes[[6]],'total watched'=outcomes[[7]])
     
# (c) dem v. rep v. pos v. neg
# outcome matrices for dem v. rep v. pos v. neg
  

outcomes=list()        
for(k in 1:7){
	outcomes[[k]]=matrix(NA,3,4)  
	colnames(outcomes[[k]])=c('Mean McAuliffe & Positive','Mean McAuliffe & Negative','Mean Cuccinelli & Positive','Mean Cuccinelli & Negative')
	rownames(outcomes[[k]])=c('Democrats','Independents','Republicans')	
}
        
dop=which(pid_lean==-1 &  dem_pos==T); don=which(pid_lean==-1 &  dem_neg==T)
iop=which(pid_lean== 0 &  dem_pos==T); ion=which(pid_lean== 0 &  dem_neg==T)
rop=which(pid_lean== 1 &  dem_pos==T); ron=which(pid_lean== 1 &  dem_neg==T)
drp=which(pid_lean==-1 & rep_pos==T); drn=which(pid_lean==-1 & rep_neg==T)
irp=which(pid_lean== 0 & rep_pos==T); irn=which(pid_lean== 0 & rep_neg==T)
rrp=which(pid_lean== 1 & rep_pos==T); rrn=which(pid_lean== 1 & rep_neg==T)
 
for(k in 1:7){
	outcomes[[k]][,1] = c(mean(ymat[dop,k],na.rm=T),mean(ymat[iop,k],na.rm=T),mean(ymat[rop,k],na.rm=T))
	outcomes[[k]][,2] = c(mean(ymat[don,k],na.rm=T),mean(ymat[ion,k],na.rm=T),mean(ymat[ron,k],na.rm=T))	   
	outcomes[[k]][,3] = c(mean(ymat[drp,k],na.rm=T),mean(ymat[irp,k],na.rm=T),mean(ymat[rrp,k],na.rm=T))
	outcomes[[k]][,4] = c(mean(ymat[drn,k],na.rm=T),mean(ymat[irn,k],na.rm=T),mean(ymat[rrn,k],na.rm=T))	
}

cand_tone_by_pid=list('skipped'=outcomes[[1]],'replay'=outcomes[[2]],'share'=outcomes[[3]],'getlink'=outcomes[[4]],'all'=outcomes[[5]],'time watched'=outcomes[[6]],'total watched'=outcomes[[7]])



# PID STRATIFICATION; like randomizing within blocks

set.seed(1005)

dems=inds=reps=cand_dems=cand_inds=cand_reps=tone_dems=tone_inds=tone_reps=tone_partisan=cand_partisan=list()               
for(k in 1:7){
	dems[[k]]=inds[[k]]=reps[[k]]=matrix(NA,2000,4)
	tone_dems[[k]]=tone_inds[[k]]=tone_reps[[k]]=tone_partisan[[k]]=cand_partisan[[k]]=cand_dems[[k]]=cand_inds[[k]]=cand_reps[[k]]=matrix(NA,2000,2)  
	colnames(dems[[k]])=colnames(inds[[k]])=colnames(reps[[k]])=c('Mean McAuliffe & Positive','Mean McAuliffe & Negative','Mean Cuccinelli & Positive','Mean Cuccinelli & Negative')
	colnames(cand_partisan[[k]])=colnames(cand_dems[[k]])=colnames(cand_inds[[k]])=colnames(cand_reps[[k]])=c('Mean McAuliffe','Mean Cuccinelli')
	colnames(tone_partisan[[k]])=colnames(tone_dems[[k]])=colnames(tone_inds[[k]])=colnames(tone_reps[[k]])=c('Mean Positive','Mean Negative')     
}

#all_pid=skip_time_pid=skip_share_pid=skip_getlink_pid=skip_replay_pid=skip_means_pid=matrix(NA,2000,4*3)
#pvn_all=pvn_skip=pvn_replay=pvn_share=pvn_getlink=matrix(NA,2000,6)        
tr_un=sort(unique(tr))
        
#vec=array(NA,2000)
for(j in 1:2000){
	indsa=sample(1:length(tr),size=length(tr),replace=T)
	omega=tr[indsa]
	pidl=pid_lean[indsa]

	ym=ymat[indsa,]

	# pid disag
	
	ind1 = which(str_sub(omega,2,2) == "P" & str_sub(omega,1,1) == "M" & pidl==-1)
	ind2 = which(str_sub(omega,2,2) == "N" & str_sub(omega,1,1) == "M" & pidl==-1)
	ind3 = which(str_sub(omega,2,2) == "P" & str_sub(omega,1,1) == "C" & pidl==-1)
	ind4 = which(str_sub(omega,2,2) == "N" & str_sub(omega,1,1) == "C" & pidl==-1)   
	
	ind5 = which(str_sub(omega,2,2) == "P" & str_sub(omega,1,1) == "M" & pidl==0)
	ind6 = which(str_sub(omega,2,2) == "N" & str_sub(omega,1,1) == "M" & pidl==0)
	ind7 = which(str_sub(omega,2,2) == "P" & str_sub(omega,1,1) == "C" & pidl==0)
	ind8 = which(str_sub(omega,2,2) == "N" & str_sub(omega,1,1) == "C" & pidl==0)

	ind9 = which(str_sub(omega,2,2) == "P" & str_sub(omega,1,1) == "M" & pidl==1)
   ind10 = which(str_sub(omega,2,2) == "N" & str_sub(omega,1,1) == "M" & pidl==1)
   ind11 = which(str_sub(omega,2,2) == "P" & str_sub(omega,1,1) == "C" & pidl==1)
   ind12 = which(str_sub(omega,2,2) == "N" & str_sub(omega,1,1) == "C" & pidl==1)    

	# pooling over party of candidate 
	indtone1 = which(str_sub(omega,2,2) == "P" & pidl==-1)
	indtone2 = which(str_sub(omega,2,2) == "N" & pidl==-1)                                              
	indtone3 = which(str_sub(omega,2,2) == "P" & pidl== 0)
	indtone4 = which(str_sub(omega,2,2) == "N" & pidl== 0)
	indtone5 = which(str_sub(omega,2,2) == "P" & pidl== 1)
	indtone6 = which(str_sub(omega,2,2) == "N" & pidl== 1)     
	indtone7 = which(str_sub(omega,2,2) == 'P' & pidl!= 0)
	indtone8 = which(str_sub(omega,2,2) == 'N' & pidl!= 0)                                                        
    
	# pooling over tone
	indcand1 = which(str_sub(omega,1,1) == "M" & pidl==-1)
	indcand2 = which(str_sub(omega,1,1) == "C" & pidl==-1)                                              
	indcand3 = which(str_sub(omega,1,1) == "M" & pidl== 0)
	indcand4 = which(str_sub(omega,1,1) == "C" & pidl== 0)
	indcand5 = which(str_sub(omega,1,1) == "M" & pidl== 1)
	indcand6 = which(str_sub(omega,1,1) == "C" & pidl== 1)    
	indcand7 = which(str_sub(omega,1,1) == 'M' & pidl!= 0)
	indcand8 = which(str_sub(omega,1,1) == 'C' & pidl!= 0)                                                          

	for(k in 1:7){
		dems[[k]][j,1] = mean(ym[ind1,k],na.rm=T)
		dems[[k]][j,2] = mean(ym[ind2,k],na.rm=T)
		dems[[k]][j,3] = mean(ym[ind3,k],na.rm=T)
		dems[[k]][j,4] = mean(ym[ind4,k],na.rm=T) 
		
		inds[[k]][j,1] = mean(ym[ind5,k],na.rm=T)
		inds[[k]][j,2] = mean(ym[ind6,k],na.rm=T)
		inds[[k]][j,3] = mean(ym[ind7,k],na.rm=T)
		inds[[k]][j,4] = mean(ym[ind8,k],na.rm=T)
		
		reps[[k]][j,1] = mean(ym[ind9,k],na.rm=T)
	    reps[[k]][j,2] = mean(ym[ind10,k],na.rm=T)
	    reps[[k]][j,3] = mean(ym[ind11,k],na.rm=T)
	    reps[[k]][j,4] = mean(ym[ind12,k],na.rm=T) 
		
	   	tone_dems[[k]][j,1] = mean(ym[indtone1,k],na.rm=T)
		tone_dems[[k]][j,2] = mean(ym[indtone2,k],na.rm=T)

		tone_inds[[k]][j,1] = mean(ym[indtone3,k],na.rm=T)
		tone_inds[[k]][j,2] = mean(ym[indtone4,k],na.rm=T)

	    tone_reps[[k]][j,1] = mean(ym[indtone5,k],na.rm=T)
		tone_reps[[k]][j,2] = mean(ym[indtone6,k],na.rm=T)   
		        
		tone_partisan[[k]][j,1] = mean(ym[indtone7,k],na.rm=T)
		tone_partisan[[k]][j,2] = mean(ym[indtone8,k],na.rm=T)
		
		cand_dems[[k]][j,1] = mean(ym[indcand1,k],na.rm=T)
		cand_dems[[k]][j,2] = mean(ym[indcand2,k],na.rm=T)

		cand_inds[[k]][j,1] = mean(ym[indcand3,k],na.rm=T)
		cand_inds[[k]][j,2] = mean(ym[indcand4,k],na.rm=T)

	    cand_reps[[k]][j,1] = mean(ym[indcand5,k],na.rm=T)
		cand_reps[[k]][j,2] = mean(ym[indcand6,k],na.rm=T)
		
		cand_partisan[[k]][j,1] = mean(ym[indcand7,k],na.rm=T)
		cand_partisan[[k]][j,2] = mean(ym[indcand8,k],na.rm=T)
				
	}            		
}     


skip_means_pid=cbind(dems[[1]],inds[[1]],reps[[1]])
skip_replay_pid  =cbind(dems[[2]],inds[[2]],reps[[2]])    
skip_share_pid   =cbind(dems[[3]],inds[[3]],reps[[3]])    
skip_getlink_pid =cbind(dems[[4]],inds[[4]],reps[[4]])    
skip_all_pid     =cbind(dems[[5]],inds[[5]],reps[[5]])    
skip_time_pid    =cbind(dems[[6]],inds[[6]],reps[[6]])    
skip_total_pid   =cbind(dems[[7]],inds[[7]],reps[[7]])    
  
pvn_skip 	=cbind(tone_dems[[1]],tone_inds[[1]],tone_reps[[1]])    
pvn_replay	=cbind(tone_dems[[2]],tone_inds[[2]],tone_reps[[2]])    
pvn_share 	=cbind(tone_dems[[3]],tone_inds[[3]],tone_reps[[3]])    
pvn_getlink =cbind(tone_dems[[4]],tone_inds[[4]],tone_reps[[4]])    
pvn_all    	=cbind(tone_dems[[5]],tone_inds[[5]],tone_reps[[5]])    
pvn_time  	=cbind(tone_dems[[6]],tone_inds[[6]],tone_reps[[6]])    
pvn_total 	=cbind(tone_dems[[7]],tone_inds[[7]],tone_reps[[7]])    
                                                                  
pvn_skip_alt 	=cbind(tone_partisan[[1]],tone_inds[[1]])    
pvn_replay_alt	=cbind(tone_partisan[[2]],tone_inds[[2]])
pvn_share_alt 	=cbind(tone_partisan[[3]],tone_inds[[3]])
pvn_getlink_alt =cbind(tone_partisan[[4]],tone_inds[[4]])
pvn_all_alt    	=cbind(tone_partisan[[5]],tone_inds[[5]])
pvn_time_alt  	=cbind(tone_partisan[[6]],tone_inds[[6]])
pvn_total_alt 	=cbind(tone_partisan[[7]],tone_inds[[7]])


ovr_skip 	=cbind(cand_dems[[1]],cand_inds[[1]],cand_reps[[1]])    
ovr_replay  =cbind(cand_dems[[2]],cand_inds[[2]],cand_reps[[2]])    
ovr_share   =cbind(cand_dems[[3]],cand_inds[[3]],cand_reps[[3]])    
ovr_getlink =cbind(cand_dems[[4]],cand_inds[[4]],cand_reps[[4]])    
ovr_all     =cbind(cand_dems[[5]],cand_inds[[5]],cand_reps[[5]])    
ovr_time    =cbind(cand_dems[[6]],cand_inds[[6]],cand_reps[[6]])    
ovr_total   =cbind(cand_dems[[7]],cand_inds[[7]],cand_reps[[7]])    
  
ovr_skip_alt 	=cbind(cand_partisan[[1]],cand_inds[[1]])
ovr_replay_alt  =cbind(cand_partisan[[2]],cand_inds[[2]])
ovr_share_alt   =cbind(cand_partisan[[3]],cand_inds[[3]])
ovr_getlink_alt =cbind(cand_partisan[[4]],cand_inds[[4]])
ovr_all_alt     =cbind(cand_partisan[[5]],cand_inds[[5]])
ovr_time_alt    =cbind(cand_partisan[[6]],cand_inds[[6]])
ovr_total_alt   =cbind(cand_partisan[[7]],cand_inds[[7]])

skip_title='Rate of Skipping'
replay_title='Rate of Replaying'
share_title='Rate of Sharing'
getlink_title='Rate of Getting Linked'
all_title='Total Ad-Seeking Choices'              
time_title='Time Spent Watching Ad'
total_title='Total Time Watching'       

{
# summary centipede plot of four seeking behaviors
est=c(colMeans(1-skip_means_pid[,c(1:4,9:12)]),colMeans(skip_replay_pid[,c(1:4,9:12)]),colMeans(skip_share_pid[,c(1:4,9:12)]),colMeans(skip_getlink_pid[,c(1:4,9:12)]))
ci5=c(apply(1-skip_means_pid[,c(1:4,9:12)],2,quantile,prob=.025),apply(skip_replay_pid[,c(1:4,9:12)],2,quantile,prob=.025),apply(skip_share_pid[,c(1:4,9:12)],2,quantile,prob=.025),apply(skip_getlink_pid[,c(1:4,9:12)],2,quantile,prob=.025))
ci95=c(apply(1-skip_means_pid[,c(1:4,9:12)],2,quantile,prob=.975),apply(skip_replay_pid[,c(1:4,9:12)],2,quantile,prob=.975),apply(skip_share_pid[,c(1:4,9:12)],2,quantile,prob=.975),apply(skip_getlink_pid[,c(1:4,9:12)],2,quantile,prob=.975))
ns=t(c(table(interaction(tr_full[which(pid_lean!=0)],pid_lean[which(pid_lean!=0)]))))   
ns=c(ns,ns,ns,ns)
segs=rbind(est,ci5,ci95,ns)

rownames(segs) <- c("Estimate", "Quantile .025", "Quantile .975", "N")
colnames(segs) <- c( 
	'McAuliffe & Positive ','McAuliffe & Negative ','Cuccinelli & Positive ','Cuccinelli & Negative ',
	 'McAuliffe & Positive ','McAuliffe & Negative ','Cuccinelli & Positive ','Cuccinelli & Negative ',	
	'McAuliffe & Positive ','McAuliffe & Negative ','Cuccinelli & Positive ','Cuccinelli & Negative ',
	 'McAuliffe & Positive ','McAuliffe & Negative ','Cuccinelli & Positive ','Cuccinelli & Negative ',                                                                                                                                    
	'McAuliffe & Positive ','McAuliffe & Negative ','Cuccinelli & Positive ','Cuccinelli & Negative ',
	 'McAuliffe & Positive ','McAuliffe & Negative ','Cuccinelli & Positive ','Cuccinelli & Negative ',
	'McAuliffe & Positive ','McAuliffe & Negative ','Cuccinelli & Positive ','Cuccinelli & Negative ',
	 'McAuliffe & Positive ','McAuliffe & Negative ','Cuccinelli & Positive ','Cuccinelli & Negative '

)

segs=segs[,dim(segs)[2]:1]     

centipede.plot(segs, sort.segs=FALSE, main="Mean Rates of Ad Seeking Behaviors", 
	xlab="Bootstrapped Means and 95 Percent Confidence Intervals", vgrid=0, mar=c(5,12,3,5), lty=1,lwd=2, 
	right.labels=paste(' ',round(segs[1, ], 2), " (", segs[4, ], ")", sep = ""), bg=c(rep('red',4),rep('darkblue',4)), ylim=c(1,33),xlim=c(-.03,.77), col=c(rep('red',4),rep('navyblue',4)))
abline(v=50, lty=2, col="grey")  
abline(h=32.5, lty=1, col="grey")  
#abline(h=c(4.5, 8.5, 12.5, 16.5,20.5,24.5,28.5), col="grey",lty=2)
abline(h=c(8.5, 16.5,24.5), col="grey",lty=2)
text( c(-.03,-.03,-.03,-.03) ,c(4.5,12.5,20.5,28.5), labels=c("Get Link", "Share", "Replay", "Not Skip"), srt=90,cex=1.2)
legend(x=.03,y=35,
 legend=c("Democrats"),
 col=c("navyblue"), 
 lwd=3.15, 
 bty="n", 
 cex=1.2, 
 inset=.00)      

legend(x=.37,y=35,
 legend=c("Republicans"),
 col=c("red"), 
 lwd=3.15, 
 bty="n", 
 cex=1.2, 
 inset=.00)

dev.print(pdf, file="~/Dropbox/Seeing_Spots/replication/figures/GovSummaryInteractions_60.pdf", width=11, height=13,pointsize=12)
}

# END figure_xii.R